!!!! Business_Cycles_Insurance Toolbox - (12/21/2022) !!!!
!! Author: Michael Irwin
!! Subroutines for "The Impact of Private and Public Insurance on Consumption Over the Business Cycle"

!!!! Table of Contents !!!!
!! 1.) Forecasting Rules
!! 2.) Aggregate Fluctuations
!! 3.) UI Polices
!! 4.) Employment Transition Matrix
!! 5.) Individual State Transition Matrix
!! 6.) Aggregate State Transition Matrix
!! 7.) Utility Function (CRRA)
!! 8.) Borrowing Limit
!! 9.) Grid Search Subroutine

module BCI_TB

use BCI_Par
use Fortran_TB

implicit none

    contains

!!Introduction to Code
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 1.) Forecasting Rules 
    
subroutine forecasting_rules(switch_scenario, avg_K, sigma_K, nu_K0, nu_K1, nu_K2, avg_L, sigma_L, nu_L0, nu_L1, nu_L2)

    !!!! Variables and Arrays !!!!
    integer, intent(in) :: switch_scenario
    real(rp), intent(out) :: avg_K, sigma_K, avg_L, sigma_L
    real(rp), dimension(Nz), intent(out) :: nu_K0, nu_K1, nu_K2, nu_L0, nu_L1, nu_L2

    !!!! Forecasting Rules - Use Benchmark Forecasting as initial guess !!!!
    avg_K = 37.87827404_rp
    sigma_K = 0.37458615_rp
    nu_k0(1) = 0.03663775_rp
    nu_k1(1) = 0.98449003_rp
    nu_k2(1) = 0.07347824_rp
    nu_k0(2) = 0.03353407_rp
    nu_k1(2) = 0.98527738_rp
    nu_k2(2) = 0.07649407_rp
    
    !avg_K = 37.79133805_rp
    !sigma_K = 0.39779416_rp
    !nu_k0(1) = 0.06024371_rp
    !nu_k1(1) = 0.97863030_rp
    !nu_k2(1) = 0.06550780_rp
    !nu_k0(2) = 0.05570280_rp
    !nu_k1(2) = 0.98001571_rp
    !nu_k2(2) = 0.06341740_rp
    
    avg_L = 1.30661970_rp
    sigma_L = 0.01594179_rp
    nu_l0(1) = 0.38151479_rp
    nu_l1(1) = 0.32025702_rp
    nu_l2(1) = -0.05384542_rp
    nu_l0(2) = 0.24206183_rp
    nu_l1(2) = 0.48992780_rp
    nu_l2(2) = -0.03363227_rp
    
    !avg_L = 1.30387056_rp
    !sigma_L = 0.01623496_rp
    !nu_l0(1) = 0.25150517_rp
    !nu_l1(1) = 0.32582389_rp
    !nu_l2(1) = -0.01887417_rp
    !nu_l0(2) = 0.14849854_rp
    !nu_l1(2) = 0.50046075_rp
    !nu_l2(2) = -0.00896964_rp


end subroutine forecasting_rules
    
!! 1.) Forecasting Rules
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 2.) Aggregate Fluctuations 
    
subroutine aggregate_fluctuations(switch_scenario, zgrid, lambda, xi, iota)

!!!! Variables and Arrays !!!!
integer, intent(in) :: switch_scenario
real(rp), dimension(Nz), intent(out) :: zgrid, lambda, xi, iota

!!!! Aggregate Fluctuations !!!!
zgrid(1) = z_g
lambda(1) = lambda_g
xi(1) = xi_g
iota(1) = iota_g
iota(2) = iota_b
if ( switch_scenario == 1 ) then
    zgrid(2) = z_b
    lambda(2) = lambda_b
    xi(2) = xi_b
elseif ( switch_scenario == 2 ) then            !! TFP Decomposition
    zgrid(2) = z_g
    lambda(2) = lambda_b
    xi(2) = xi_b
elseif ( switch_scenario == 3 ) then            !! Separation Decomposition
    zgrid(2) = z_b
    lambda(2) = lambda_b
    xi(2) = xi_g
elseif ( switch_scenario == 4 ) then            !! Finding Decomposition
    zgrid(2) = z_b
    lambda(2) = lambda_g
    xi(2) = xi_b
elseif (switch_scenario >= 5 .and. switch_scenario <= 12 ) then
    zgrid(2) = z_b
    lambda(2) = lambda_b
    xi(2) = xi_b
elseif ( switch_scenario == 13 ) then           !! Recalibrate: No job finding rate fluctuations
    zgrid(2) = z_b
    lambda(2) = lambda_g
    xi(2) = 0.064_rp
end if


end subroutine aggregate_fluctuations
    
!! 2.) Aggregate Fluctuations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 3.) UI Policies

subroutine UI_policies(switch_scenario, upsilon_d, upsilon_r, upsilon_bar)

!!!! Variables and Arrays !!!!
integer, intent(in) :: switch_scenario
real(rp), dimension(Nz), intent(out) :: upsilon_d, upsilon_r, upsilon_bar

!!!! UI Policies !!!!
upsilon_d(1) = upsilon_dg
upsilon_r(1) = upsilon_rg
upsilon_bar(1) = upsilon_bar_g
if ( switch_scenario <= 4 .or. switch_scenario == 7 .or. switch_scenario == 10 .or. switch_scenario == 13 ) then        !! UI Extension
    upsilon_d(2) = upsilon_db
    upsilon_r(2) = upsilon_rg
    upsilon_bar(2) = upsilon_bar_g
elseif ( switch_scenario == 5 .or. switch_scenario == 8 .or. switch_scenario == 11 ) then                               !! Acyclical UI
    upsilon_d(2) = upsilon_dg
    upsilon_r(2) = upsilon_rg
    upsilon_bar(2) = upsilon_bar_g
elseif (switch_scenario == 6 .or. switch_scenario == 9 .or. switch_scenario == 12 ) then                                !! UI Level Increase
    upsilon_d(2) = upsilon_dg
    upsilon_r(2) = upsilon_rb
    upsilon_bar(2) = upsilon_bar_b
end if

end subroutine UI_policies

!! 3.) UI Policies
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 4.) Employment Transition Matrix

subroutine emp_transitions(n, xi, lambda, upsilon_d, pi)

!!!! Variables and Arrays !!!!
integer, intent(in) :: n
real(rp), intent(in) :: xi, lambda, upsilon_d
real(rp), dimension(n,n), intent(out) :: pi

!!!! Transition Matrix !!!!
pi = zero
pi(1,1) = one - xi
pi(1,2) = xi
pi(1,3) = zero
pi(2,1) = lambda
pi(2,2) = (one - upsilon_d)*(one - lambda)
pi(2,3) = upsilon_d*(one - lambda)
pi(3,1) = lambda
pi(3,2) = zero
pi(3,3) = one - lambda

end subroutine emp_transitions

!! 4.) Employment Transition Matrix
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 5.) Individual Transition Matrix

subroutine ind_transitions(ni, ne, nn, nb, pi_e, pi_n, e_index, n_index, b_index, i_index, pi_i )

!!!! Variables and Arrays !!!!
integer, intent(in) :: ni, ne, nn, nb
real(rp), dimension(ne,ne), intent(in) :: pi_e
real(rp), dimension(nn,nn), intent(in) :: pi_n
integer :: i, e, n, b, i2, e2, n2, b2, count
integer, dimension(ni), intent(out) :: e_index, n_index, b_index
integer, dimension(ne,nn,nb), intent(out) :: i_index
real(rp), dimension(ni,ni,2), intent(out) :: pi_i

!!!! Assign Indices !!!!
count = 0
do e = 1,ne
	do n = 1,nn
		do b = 1,nb
			count = count + 1
			e_index(count) = e
			n_index(count) = n
			b_index(count) = b
            i_index(e,n,b) = count
		end do
	end do
end do

!!!! Overall Transition Matrix !!!!
pi_i = zero
do i = 1,ni
	n = n_index(i)
	e = e_index(i)
	b = b_index(i)

	do i2 = 1,ni
		n2 = n_index(i2)
		e2 = e_index(i2)
		b2 = b_index(i2)
        
		if (b == b2)  then
			if ( n == 1 ) then
				if ( n2 == 1 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)*pi_e(e,e2)
                elseif ( n2 == 2 .and. e == e2 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)
                elseif (n2 == 3 .and. e == e2 ) then
                    pi_i(i,i2,2) = one
                endif
            elseif (n == 2 ) then
                if ( n2 == 1 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)*pi_e(e,e2)
                elseif ( n2 == 2 .and. e == e2) then
                    pi_i(i,i2,1) = pi_n(n,n2)
                elseif ( n2 == 3 .and. e == e2 )  then
                    pi_i(i,i2,1) = pi_n(n,n2)
                    pi_i(i,i2,2) = one
                end if
            else
                if ( n2 == 1 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)*pi_e(e,e2)
                elseif ( n2 == 3 .and. e == e2 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)
                    pi_i(i,i2,2) = one
                end if
            end if
		end if
        
	end do
end do

end subroutine ind_transitions

!! 5.) Individual Transition Matrix
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 6.) Aggregate Transition Matrix

subroutine agg_transitions( nz, rho_z1, rho_z2, pi_z, cum_z )

!!!! Variables and Arrays !!!!
integer, intent(in) :: nz
real(rp), intent(in) :: rho_z1, rho_z2
real(rp), dimension(nz,nz), intent(out) :: pi_z, cum_z

!!!! Aggregate Transition Matrix !!!!
pi_z = zero
cum_z = zero
pi_z(1,1) = rho_z1
pi_z(1,2) = one - rho_z1
pi_z(2,1) = one - rho_z2
pi_z(2,2) = rho_z2
cum_z(1,1) = pi_z(1,1)
cum_z(1,2) = pi_z(1,1) + pi_z(1,2)
cum_z(2,1) = pi_z(2,1)
cum_z(2,2) = pi_z(2,1) + pi_z(2,2)

end subroutine agg_transitions

!! 6.) Aggregate Transitions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 7.) Utility Function

subroutine utility(beta, chi, consumption, EV, u_val)

!!!! Variables and Arrays !!!!
real(rp), intent(in) :: beta, chi, consumption, EV
real(rp), intent(out) :: u_val

!!!! Function !!!!
u_val = (one/(one-sigma))*(( consumption )**(one-sigma)) - chi + beta*EV

end subroutine utility

!! 7.) Utility Function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 8.) Borrowing Limit

subroutine borrowing_limit(nd, xgrid, discount, iB, minB)

!!!! Variables and Arrays !!!!
integer, intent(in) :: nd
real(rp), dimension(nd), intent(in) :: xgrid, discount
integer, intent(out) :: iB
real(rp), intent(out) :: minB
integer :: x

!!!! Calculate Endogenous Borrowing Limit !!!!
minB = 0.0_rp
iB = nd
do x = 1,nd
	if ( xgrid(x)*discount(x) < minB ) then
		minB = xgrid(x)*discount(x)
		iB = x
	end if
end do

end subroutine borrowing_limit

!! 8.) Borrowing Limit
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!! 9.) Grid Search

subroutine grid_search(ibar, beta, chi, a_val, income, discount, agrid, EV, maxV, maxG, maxC, maxQ)

!!!! Variables and Arrays !!!!
integer, intent(in) :: ibar
real(rp), intent(in) :: beta, chi, a_val, income
real(rp), dimension(nd), intent(in) :: discount
real(rp), dimension(na), intent(in) :: agrid, EV
real(rp), intent(out) :: maxV, maxG, maxC, maxQ
integer :: a2
real(rp) :: consumption, Value

!!!! Optimal Value for Borrowing !!!!
maxV = disutility; maxG = zero; maxC = zero; maxQ = zero
do a2 = ibar,nd
    consumption = a_val + income + y_bar - discount(a2)*agrid(a2)
    if ( consumption > zero ) then
        Value = (one/(one-sigma))*(( consumption )**(one-sigma)) - chi + beta*EV(a2)
    else
        Value = disutility
    end if

    if ( Value > maxV ) then
    	maxV = Value
    	maxG = agrid(a2)
    	maxC = consumption
        maxQ = discount(a2)
    end if
end do

end subroutine grid_search

!! 9.) Grid Search
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!! 10.) Income Distribution

subroutine income_distribution(nx, e_index, n_index, upsilon_r, gamma_j, egrid, i_income, j_income, income_grid)

!!!! Variables and Arrays !!!!
integer :: i, i2, j, e, n, count, ix
real(rp) :: income, omega_x
integer, intent(in) :: nx       !!dimension of income grid (ni*nj)
integer, dimension(ni), intent(in) :: e_index, n_index
real(rp), intent(in) :: upsilon_r
real(rp), dimension(nj), intent(in) :: gamma_j
real(rp), dimension(ne), intent(in) :: egrid
integer, dimension(nx), intent(out) :: i_income, j_income
real(rp), dimension(nx) :: income_grid

!!!! Initialize Income Grid !!!!
i_income(:) = 1
j_income(:) = 1
income_grid(:) = zero

!!!! Income Order - Positive Income !!!!
count = 1
do j = 1,nj
    do i = 1,ni
        e = e_index(i); n = n_index(i)
        
        if ( n == 1 ) then
            income = gamma_j(j)*egrid(e)
        elseif ( n == 2 ) then
            income = upsilon_r*gamma_j(j)*egrid(e)
        else
            income = zero
        end if
        
        if ( count < 3 ) then
            i_income(count) = i
            j_income(count) = j
            income_grid(count) = income
            count = count + 1
        else
            
            if ( income >= income_grid(count-1) ) then
                i_income(count) = i
                j_income(count) = j
                income_grid(count) = income
                count = count + 1
            elseif ( income <= income_grid(1) ) then
                do i2 = count,2,-1
                    i_income(i2) = i_income(i2-1)
                    j_income(i2) = j_income(i2-1)
                    income_grid(i2) = income_grid(i2-1)
                end do
                i_income(1) = i
                j_income(1) = j
                income_grid(1) = income
                count = count + 1
            else
                call interpolation(income_grid, income, 1, count-1, nx, ix, omega_x)
                do i2 = count,ix+2,-1
                    income_grid(i2) = income_grid(i2-1)
                    i_income(i2) = i_income(i2-1)
                    j_income(i2) = j_income(i2-1)
                end do
                income_grid(ix+1) = income
                i_income(ix+1) = i
                j_income(ix+1) = j
                count = count + 1
            end if
            
        end if
        
    end do
end do


end subroutine income_distribution

!! 10.) Income Distribution
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!! 11.) Future Index - For Simulation

subroutine future_index(cum_pi, rand_x, nx, ix)

!!!! Variables and Arrays !!!!
integer :: x
real(rp) :: omega
integer, intent(in) :: nx
real(rp), intent(in) :: rand_x
real(rp), dimension(nx), intent(in) :: cum_pi
integer, intent(out) :: ix

!!!! Calculate Future Index !!!!
if ( rand_x <= cum_pi(1) ) then
    ix = 1
else
    call interpolation(cum_pi, rand_x, 1, nx, nx, x, omega)
    ix = x + 1
end if

end subroutine


end module BCI_TB